﻿Imports System.Collections.Generic
Imports System.Math

''' <summary>
''' Routines to convert between UTM and Lat/Long coordinates
''' </summary>
''' <remarks>Based on spreadsheets obtained from Ivan Dominguez, ACP, via www.geofumadas.com</remarks>
Public Class clsCoordinates
    Private Class clsEllipsoid
        Friend Name As String
        Friend A As Double 'major axis
        Friend B As Double 'minor axis
    End Class

    Public Enum enumEllipsoid
        Airy_1830
        Airy_Modificado_1965
        Bessel_1841
        Clarke_1866
        Clarke_1880
        Fischer_1968
        GRS80
        Hayford_1909
        Helmert_1906
        Hough_1960
        Internacional_1909
        Internacional_1924
        Krasovsky_1940
        Mercury_1960
        Mercury_Modificado_1968
        Nuevo_International_1967
        Sudamericano_1969
        Walbeck_1817
        WGS66
        WGS72
        WGS84
    End Enum

    Public Shared Function Ellipsoids() As List(Of String)
        Dim lst As New List(Of String)
        For Each s As String In [Enum].GetNames(GetType(enumEllipsoid))
            lst.Add(s.Replace("_", " "))
        Next
        Return lst
    End Function

    Public Enum enumNorthSouth
        North = 1
        South = -1
    End Enum

    Public Enum enumEastWest
        East = 1
        West = -1
    End Enum

    Private Ellipsoid As enumEllipsoid
    Private Hemisphere As enumNorthSouth
    Private Zone As Integer

    Private dictEllipsoid As New Dictionary(Of enumEllipsoid, clsEllipsoid)

    Private f, n, AA, a(3), b(3), d(3), N0, E0, k0 As Double

    Public Sub New(Ellipsoid As enumEllipsoid, Zone As Integer, Hemisphere As enumNorthSouth)
        With dictEllipsoid
            .Add(enumEllipsoid.Airy_1830, New clsEllipsoid With {.Name = "Airy 1830", .A = 6377563.396, .B = 6356256.91})
            .Add(enumEllipsoid.Airy_Modificado_1965, New clsEllipsoid With {.Name = "Airy Modificado 1965", .A = 6377340.189, .B = 6356034.4479})
            .Add(enumEllipsoid.Bessel_1841, New clsEllipsoid With {.Name = "Bessel 1841", .A = 6377397.155, .B = 6356078.96284})
            .Add(enumEllipsoid.Clarke_1866, New clsEllipsoid With {.Name = "Clarke 1866", .A = 6378206.4, .B = 6356583.8})
            .Add(enumEllipsoid.Clarke_1880, New clsEllipsoid With {.Name = "Clarke 1880", .A = 6378249.145, .B = 6356514.86955})
            .Add(enumEllipsoid.Fischer_1968, New clsEllipsoid With {.Name = "Fischer 1968", .A = 6378150.0, .B = 6356768.33})
            .Add(enumEllipsoid.GRS80, New clsEllipsoid With {.Name = "GRS80", .A = 6378137.0, .B = 6356752.31414})
            .Add(enumEllipsoid.Hayford_1909, New clsEllipsoid With {.Name = "Hayford 1909", .A = 6378388.0, .B = 6356911.94613})
            .Add(enumEllipsoid.Helmert_1906, New clsEllipsoid With {.Name = "Helmert 1906", .A = 6378200.0, .B = 6356818.17})
            .Add(enumEllipsoid.Hough_1960, New clsEllipsoid With {.Name = "Hough 1960", .A = 6378270.0, .B = 6356794.343479})
            .Add(enumEllipsoid.Internacional_1909, New clsEllipsoid With {.Name = "Internacional 1909", .A = 6378388.0, .B = 6356911.94613})
            .Add(enumEllipsoid.Internacional_1924, New clsEllipsoid With {.Name = "Internacional 1924", .A = 6378388.0, .B = 6356911.94613})
            .Add(enumEllipsoid.Krasovsky_1940, New clsEllipsoid With {.Name = "Krasovsky 1940", .A = 6378245.0, .B = 6356863.0188})
            .Add(enumEllipsoid.Mercury_1960, New clsEllipsoid With {.Name = "Mercury 1960", .A = 6378166.0, .B = 6356784.283666})
            .Add(enumEllipsoid.Mercury_Modificado_1968, New clsEllipsoid With {.Name = "Mercury Modificado 1968", .A = 6378150.0, .B = 6356768.337303})
            .Add(enumEllipsoid.Nuevo_International_1967, New clsEllipsoid With {.Name = "Nuevo International 1967", .A = 6378157.5, .B = 6356772.2})
            .Add(enumEllipsoid.Sudamericano_1969, New clsEllipsoid With {.Name = "Sudamericano 1969", .A = 6378160.0, .B = 6356774.72})
            .Add(enumEllipsoid.Walbeck_1817, New clsEllipsoid With {.Name = "Walbeck 1817", .A = 6376896.0, .B = 6355834.8467})
            .Add(enumEllipsoid.WGS66, New clsEllipsoid With {.Name = "WGS66", .A = 6378145.0, .B = 6356759.769356})
            .Add(enumEllipsoid.WGS72, New clsEllipsoid With {.Name = "WGS72", .A = 6378135.0, .B = 6356750.519915})
            .Add(enumEllipsoid.WGS84, New clsEllipsoid With {.Name = "WGS84", .A = 6378137.0, .B = 6356752.314245})
        End With
        Me.Ellipsoid = Ellipsoid
        Me.Zone = Zone
        Me.Hemisphere = Hemisphere
        If Zone < 1 Or Zone > 60 Then Throw New ArgumentOutOfRangeException("Zone", "Zone must be an integer between 1 and 60")

        With dictEllipsoid(Ellipsoid)
            'convert to km
            .A /= 1000
            .B /= 1000
            f = (.A - .B) / .A
            n = f / (2 - f)
            AA = .A / (1 + n) * (1 + n ^ 2 / 4 + n ^ 4 / 64)
            a(1) = 1.0 / 2.0 * n - 2.0 / 3.0 * n ^ 2 + 5.0 / 16.0 * n ^ 3
            a(2) = 13.0 / 48.0 * n ^ 2 - 3.0 / 5.0 * n ^ 3
            a(3) = 61.0 / 240.0 * n ^ 3
            b(1) = 1.0 / 2.0 * n - 2.0 / 3.0 * n ^ 2 + 37.0 / 96.0 * n ^ 3
            b(2) = 1.0 / 48.0 * n ^ 2 + 1.0 / 15.0 * n ^ 3
            b(3) = 17.0 / 480.0 * n ^ 3
            d(1) = 2.0 * n - 2.0 / 3.0 * n ^ 2 - 2.0 * n ^ 3
            d(2) = 7.0 / 3.0 * n ^ 2 - 8.0 / 5.0 * n ^ 3
            d(3) = 56.0 / 15.0 * n ^ 3
            N0 = IIf(Hemisphere = enumNorthSouth.North, 0, 10000)
            E0 = 500
            k0 = 0.9996
        End With
    End Sub

    ''' <summary>
    ''' Convert UTM coordinates to lat/long
    ''' </summary>
    ''' <param name="UTM_E">X coordinate in meters</param>
    ''' <param name="UTM_N">Y coordinate in meters</param>
    ''' <returns>Decimal lat/long</returns>
    ''' <remarks>See http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system </remarks>
    Public Function UTMtoLatLon(UTM_E As Double, UTM_N As Double) As PointF
        'convert from m to km
        UTM_E /= 1000
        UTM_N /= 1000
        With dictEllipsoid(Ellipsoid)
            Dim EE As Double = (UTM_N - N0) / k0 / AA
            Dim NN As Double = (UTM_E - E0) / k0 / AA
            Dim EEP As Double = EE
            Dim NNP As Double = NN
            Dim SSP As Double = 1.0
            Dim TTP As Double = 0.0
            For j As Integer = 1 To 3
                EEP -= b(j) * Sin(2.0 * j * EE) * Cosh(2.0 * j * NN)
                NNP -= b(j) * Cos(2.0 * j * EE) * Sinh(2.0 * j * NN)
                SSP -= 2.0 * j * b(j) * Cos(2.0 * j * EE) * Cosh(2.0 * j * NN)
                TTP += 2.0 * j * b(j) * Sin(2.0 * j * EE) * Sinh(2.0 * j * NN)
            Next
            Dim X As Double = Asin(Sin(EEP) / Cosh(NNP))
            Dim Lat As Double = X
            For j As Integer = 1 To 3
                Lat += d(j) * Sin(2.0 * j * X)
            Next
            Lat *= 180 / PI
            Dim Lon0 As Double = Zone * 6 - 183
            Dim Lon As Double = Lon0 + Atan(Sinh(NNP) / Cos(EEP)) * 180 / PI
            Dim k As Double = k0 * AA / .A * Sqrt((1 + ((1 - N) / (1 + N) * Tan(Lat)) ^ 2) * (Cos(EEP) ^ 2 + Sinh(NNP) ^ 2) / (SSP ^ 2 + TTP ^ 2)) 'point scale factor
            Dim c As Double = Hemisphere * Atan((TTP + SSP * Tan(EEP) * Tanh(NNP)) / (SSP - TTP * Tan(EEP) * Tanh(NNP))) 'meridian convergence
            Return New PointF(Lon, Lat)
        End With
    End Function

    Private Function Atanh(x As Double) As Double
        'see http://msdn.microsoft.com/en-us/library/w3t84e33.aspx
        Return Log((1 + x) / (1 - x)) / 2
    End Function

    ''' <summary>
    ''' Convert latitude and longitude to UTM coordinates
    ''' </summary>
    ''' <param name="Lon">Longitude; if negative will assume is west and will ignore EW parameter</param>
    ''' <param name="Lat">Latitude; N/S is taken from the UTM hemisphere</param>
    ''' <param name="EW">If Longitude is positive, can set this to West to force west longitude</param>
    ''' <returns>X,Y coordinates in UTM zone (m)</returns>
    Public Function LatLontoUTM(Lon As Double, Lat As Double, Optional EW As enumEastWest = enumEastWest.East) As PointF
        With dictEllipsoid(Ellipsoid)
            If EW = enumEastWest.West Then Lon = -Abs(Lon) 'force it to be negative longitude
            Dim Lon0 As Double = Zone * 6 - 183
            Dim t As Double = Sinh(Atanh(Sin(Lat) - 2 * Sqrt(n) / (1 + n) * Atanh(2 * Sqrt(n) / (1 + n) * Sin(Lat))))
            Dim EEP As Double = Atan(t / Cos(Lon - Lon0))
            Dim NNP As Double = Atanh(Sin(Lon - Lon0) / Sqrt(1 + t ^ 2))
            Dim SS As Double = 1.0
            Dim TT As Double = 0.0
            For j As Integer = 1 To 3
                SS += 2 * j * a(j) * Cos(2 * j * EEP) * Cosh(2 * j * NNP)
                TT += 2 * j * a(j) * Sin(2 * j * EEP) * Sinh(2 * j * NNP)
            Next
            Dim esum As Double = NNP
            Dim nsum As Double = EEP
            For j As Integer = 1 To 3
                esum += a(j) * Cos(2 * j * EEP) * Sinh(2 * j * NNP)
                nsum += a(j) * Sin(2 * j * EEP) * Cosh(2 * j * NNP)
            Next
            Dim EE As Double = E0 + k0 * AA * esum
            Dim NN As Double = N0 + k0 * AA * nsum
            Return New PointF(EE * 1000, NN * 1000)
        End With
    End Function

    Public Function LatLontoUTM(LonDeg As Integer, LonMin As Integer, LonSec As Double, LatDeg As Integer, LatMin As Integer, LatSec As Double, Optional EW As enumEastWest = enumEastWest.East)
        Return LatLontoUTM(DMStoDeg(LonDeg, LonMin, LonSec), DMStoDeg(LatDeg, LatMin, LatSec), EW)
    End Function

    Private Function DMStoDeg(D As Integer, M As Integer, S As Double) As Double
        If M < 0 Or M > 60 Or S < 0.0 Or S > 60.0 Then Throw New ArgumentOutOfRangeException
        Return (Abs(D) + M / 60.0 + S / 3600.0) * (Sign(D))
    End Function

End Class
